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Abstract 

We present a practical numerical technique for calculating tunneling ionization 
rates from arbitrary 1-D potential wells in the presence of a linear external potential 
by determining the widths of the resonances in the spectral density, p(E), adiabatically 
connected to the field-free bound states. While this technique applies to more general 
external potentials, we focus on the ionization of electrons from atoms and molecules by 
DC electric fields, as this has an important and immediate impact on the understanding 
of the multiphoton ionization of molecules in strong laser fields. 

While the behavior of atoms and molecules in strong laser fields has been studied for 
almost two decades, new phenomena continue to arise and a unified picture for either does 
not yet exist [[[], |2], [| . However, simple models for tunneling ionization j|, |5|] have provided a 
crucial baseline for understanding atoms in strong fields, as deviations from this basic process 
in ion yields and electron energy spectra have pointed to new physics |l|, §, ||. Conversely, 
the lack of a general tunneling model for molecules has seriously impeded the interpretation 
of molecular data which clearly hints at even richer phenomena f||, [7] than for atoms due to 
the extra degrees of freedom present in a molecule. 

Ionization from even diatomic potential wells are more difficult to treat than atoms 
because the barrier to ionization can shift between the barrier separating the two wells 
and the outer barrier as a function of field strength and internuclear separation || (see 
Fig. 1). These situations are quite distinct and, indeed, this competition between barriers 
leads to interesting and complex behavior. Thus, it is critical for any tunneling ionization 
model to work for arbitrary potential wells so as to treat these different situations uniformly. 
Furthermore, a useful model for molecules in strong laser fields must be non-perturbative in 
both the rates and the Stark shifts. In this Letter we introduce a practical non-perturbative 
technique for computing tunneling ionization rates from an arbitrary 1-D potential well. 
After outlining the general method, we apply it to a single-well ('atomic') potential, yielding 
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excellent agreement with existing theoretical analyses. Then we apply our method to a 
double-well ('molecular') potential, like that shown in Fig. 1, to illustrate the generality of 
the method and to show new features which arise due to the double-well structure. However, 
we emphasize that the method is exactly the same for any arbitrary shaped 1-D potential 
well. 

The problem of determining ionization rates and Stark shifts non-perturbatively for an 
arbitrary potential in 1-D can be reduced to calculating the spectral density, p(E). For 
potentials with strictly bound states p(E) becomes a series of delta functions at the energy 
eigenvalues. Numerically integrating the Schrodinger equation (SE) determines these ener- 
gies by finding the energy at which the wavefunction satisfies the given boundary conditions. 
This simple approach cannot be extended to the question of tunneling ionization because, 
by definition, the electron can escape to infinity. Given this fact, the electron's state lies in 
the continuum; hence, all energies are allowed and simply integrating the SE at a particular 
energy provides little information, on its own. Rather, one must return to calculating p(E) 
in general. In the presence of a linear potential, the delta functions in p(E) for the bound 
states will evolve into Lorentzian-like resonances whose center positions reflect the Stark 
shift of the states and widths give their decay rates. As decay of the levels can only occur 
through ionization the width then directly gives the tunneling rate. 

The limitation to 1-D potential wells may seem overly restrictive and requires some 
discussion. It is well known that ionization rates of atoms and molecules depend strongly on 
field strength JO. For this reason ionization occurs overwhelmingly in the direction of the 
electric field, essentially reducing the problem to 1-D. For diatomic molecules, the ionization 
rate is much greater when the field is along the internuclear axis §| and, again, reduces the 
problem to 1-D. 

While a full discussion of p(E) is rarely presented in standard books on quantum mechan- 
ics, a general approach exists for calculating p(E), based on the Weyl-Titchmarsh-Kodaira 
(WTK) spectral theorem fL0| . Remarkably, this approach is extremely well suited to nu- 



merical implementation. Essentially, the spectral theorem states that the "resolution of the 
identity" can be expressed as 

/oo 
dE<j> i {x : E)<l> j {y,E)p ij {E) (1) 
i,j=±,* ~°° 

where 0i,2(X E) are two linearly independent eigenfunctions of the SE at energy E and Pij{E) 
is the spectral density matrix. Before discussing how to calculate and interpret pij(E), we 
must consider the question of normalization, as, in general, some of the states may lie in the 
continuum. 

The usual concept of a finite square integral clearly cannot be used for the normalization 
of continuum wavef unctions. In the WTK approach, normalization is imposed locally by 
setting the Wronskian of the eigenfunctions equal to one which also ensures that </>i and </>2 
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are linearly independent. Thus, we can choose the following initial conditions at any point 

x = c : 

0i(c) = l, #(c) = 0, 2 (c) = O, 2 (c) = l. (2) 

For convenience, we generally choose x — 0. The complete functions 0i(x) can now be found 
by numerically integrating the SE away from the origin using these initial conditions. 

To construct the spectral density matrix pij(E), we first find the linear combinations of 

0i and 2 , 

i>± (x, A) = <pi (x, A) + m± (A) 2 , (3) 

such that fl^dx \ip-\ 2 and dx \ip+\ 2 are both finite, where A = E + is and e > 0. In 
other words, when the energy has a positive imaginary part, ip± (x, A) must decay sufficiently 
rapidly as x — > ±oo to have a finite square integral on their respective half-line. This defines 
the functions m± (A). The spectral density matrix is then 
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Pij(E) = — limlm 

7T e^O 771_ — TXIa 



m-+raj r 
2 



(4) 



If the asymptotic wavefunctions of the external potential are known analytically, as is often 
the case, all that is needed to evaluate m± is the logarithmic derivative of the correctly 
decaying asymptotic wavefunctions of just the external potential, 9±(x,\): f±(x,X) = 



6'± (x, A) /6± (x, A) PHI . We can then solve for m± (A): 



m± (A) 



±<Pl 



/: 



±92, 



(5) 



Note that the overall normalization of the asymptotic wavefunctions is not important, as the 
logarithmic derivative is sufficient to determine m± (A). 

This entire procedure is very straightforward to implement numerically: starting with 
the boundary conditions in Eq. @, <\>\ and 02 are found by integrating the SE away from 
the origin with a standard Runge-Kutta technique at a particular complex energy, A. At 
the same time, m ± are evaluated from Eq. fl5|). Ideally, one would integrate out to ±oo to 
evaluate m± (A) exactly. Since this is not possible numerically, we simply evaluate m± (A) 
as a function of x far enough from the local potential to where it converges to a constant 
value. At this point, m± (A) have been determined as well as Pij(E) from Eq. 

To finish the problem, we must interpret the matrix in Eq. (^). The spectral density 
matrix, pij(E), is positive definite 0] and, thus, its eigenvalues Ai and A 2 are positive and 
represent scalar spectral density functions. In general, there are two functions because of the 
possibility of a two-fold degeneracy (corresponding, for example, to left and right moving 
solutions) on the whole line. However, in the case of a linear external potential (Fig. 1), one 
direction is inaccessible to a free particle and no degeneracy exists. In this case, either Ai or 
A2 will be on the order of e and is unimportant, leaving just one eigenvalue, resulting in a 
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simple scalar spectral density. For problems on the half-line, x > 0, with a physical boundary 
condition at x = 0, the spectral density matrix also reduces to a single spectral density 
function JTU| , |T2J. Finally, we must consider taking the limit e — > 0. Again, numerically, 
this is simple to specify: if the state can tunnel, it will have a "natural" width T giving the 
tunneling rate and an "artificial" width determined by e. One then decreases e until r 3> e. 
We have noticed that the numerical results are surprisingly robust to changes in the various 
numerical parameters. We attribute this to the local specification of the normalization and 
initial conditions for the numerical integration and the ease in determining the m± values. 

In order to test this method we considered tunneling ionization of a hydrogen atom by a 
DC electric field: 

Vr + or 

where F is the field strength and a = \[2 for a hydrogen atom. This 1-D form of the 
Coulomb potential ]13|] has proven to be an excellent model for strong field calculations 



because it preserves the long range 1/r dependence and has an infinite Rydberg-like series of 
bound states while removing the catastrophic singularity that occurs for 1-D. Furthermore, it 
has been used recently in direct calculations of the time-dependent Schrodinger equation for 
molecules 0. Fig. 2 shows the spectral function at a moderate field strength. The first two 
low lying states in the Rydberg series can be seen, while the higher states have broadened 
and merged together. We stress that the blurring of the higher states is a real physical effect, 
due to the lowering of the potential barrier by the external field, and is not a limitation of 
the method. In a weaker field more of the excited states would be seen as well-defined peaks. 
Fig. 3 shows a typical unbound wavefunction at the peak of a resonance. At this point, 
most of the wavefunction is "inside" the potential well, although it is still a continuum state. 
Fig. 4 shows the DC Stark shift as well as the first term in the second order perturbation 
theory expansion for the 1-D potential. Finally, Fig. 5 shows the tunneling ionization rate 
of the ground state as a function of field strength. In order to compare to a 3-D potential, 
we performed an angular average of the 1-D result: 

T{F) = — I dQT{F cos 6). (7) 

As seen in Fig. 5, this averaged result agrees quite well with the standard 3-D static 
atomic tunneling model of ADK [f|, However, the slight discrepancy at the higher intensity 
is intriguing. The ADK model depends on the ionization potential and the field strength, but 
does not take into account the DC Stark shift. While the shifts are small, ionization rates 
depend exponentially on ionization potential and the effect of the Stark shift on the ionization 
rate will increase with intensity. We recalculated the ADK rate using our calculation of the 
DC Stark shift, and, as seen in Fig. 5 the agreement between our angularly averaged model 
and the ADK result including the DC Stark shift is excellent. 
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The 'molecular' or double-well potential has a much richer spectrum in an external field 
than an atom, and we cannot give a comprehensive treatment, here. However, as an example, 
we show, in Fig. 6, the effect of the separation of the double-well potential on ionization 
rates for both the ground and first excited states. The potential consists of two 1-D Coulomb 
potentials (Eq. 6) separated by a distance R, and the external field, as in Fig. 1. Using 
exactly the same method of calculating the spectral function and finding the widths of the 
peaks associated with the ground and first excited states, we find that the ionization rate is 
strongly dependent on the well separation, becoming greater at larger separations. Moreover, 
the qualitative shape of the ionization rate vs. intensity curve depends on the separation 
with apparently no saturation setting in for the excited state at the larger separation. In 
fact, these are just the conditions for electron localization and enhanced ionization predicted 
in Ref. 0. These ionization rate curves can provide a framework to analyze experimental 
ion yield data from molecules. 

In summary, we have described a practical numerical technique for calculating spectral 
density functions for arbitrary 1-D potentials on the whole line from which we can extract 
level shifts and tunneling rates. This will be of great value for interpreting data from atoms 
and molecules in strong laser fields. 
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Figure Captions 

1. Molecular potential wells showing the importance of (a) the outer barrier and (b) the 
inner barrier. 

2. Spectral density function for a model atom in an external field. 

3. Potential in an external field with wavefunction on resonance. 

4. Comparison of the 2nd order Stark shift calculation with perturbation theory, in atomic 
units. 

5. Comparison of the angularly-averaged WTK tunneling rate model with a 3-D calcula- 
tion (ADK) from Ref. Q. 

6. Ionization rate of the ground and first excited state of a double- well potential for two 
different separations. 
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Figure 1: Molecular potential wells showing the importance of (a) the outer barrier and (b) 
the inner barrier. 
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Figure 2: Spectral density function for a model atom in an external field. 
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Figure 3: Potential in an external field with wavefunction on resonance. 
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Figure 4: Comparison of the 2nd order Stark shift calculation with perturbation theory, in 
atomic units. 
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Figure 5: Comparison of the angularly-averaged WTK tunneling rate model with a 3-D 
calculation (ADK) from Ref. (§. 
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Figure 6: Ionization rate of the ground and first excited state of a double-well potential for 
two different separations. 



13 



